function [xi0,xi1,xis,xip0,Xip,phi0,phi1,phis,Phip,Cs,Cp,Pi,Fp]=Solve_model(Paras)

 global  xb0  c0  Xi0 Phi0
%%%%%��������%%%%%%%%%%

     r=Paras(1);
    mu0=Paras(2);
    mu1=Paras(3);
    sigma=Paras(4);
    alpha=Paras(5);
    taue=Paras(6);
    taui=Paras(7);
    I=Paras(8);
    p=Paras(9);
    x=Paras(10);
    pi0=Paras(11);
    eta=Paras(12);
    
%     c=p*r/(1-taui);
    
  %  psi=mu1-mu0;
 %   mu=mu0+psi*pi0;
    
    Pi=pi0;
    
    %%%%%%%%%%��%%%%%%%%
    beta1=1/2-mu0/sigma^2+sqrt((1/2-mu0/sigma^2)^2+2*r/sigma^2);
    beta2=1/2-mu0/sigma^2-sqrt((1/2-mu0/sigma^2)^2+2*r/sigma^2);
       
    beta3=1/2-mu1/sigma^2+sqrt((1/2-mu1/sigma^2)^2+2*r/sigma^2);
    beta4=1/2-mu1/sigma^2-sqrt((1/2-mu1/sigma^2)^2+2*r/sigma^2);
    
    
    
    
%     %%%%�Ʋ�ˮƽ%%%%%
%     xb0=(beta2/(beta2-1))*((r-mu0)/r)*c;
%     xb1=(beta4/(beta4-1))*((r-mu1)/r)*c;
    
    
    %%%%%Ͷ��ˮƽ%%%%%%%%%%%%
     y0=[1 0.8]';
     y0=fsolve(@(P)FB(P,Paras),y0);
 
     xi0=roundn(y0(1,1),-8);
     c0=roundn(y0(2,1),-8);
     Xi0=xi0;
     xb0=(beta2/(beta2-1))*((r-mu0)/r)*c0;
     
     %%%%%%%������%%%%%%%%%%%%%%
     E0=(1-taue)*(xi0/(r-mu0)-c0/r+(c0/r-xb0/(r-mu0))*(xi0/xb0)^beta2);
     D0=(1-taui)*c0/r+((1-alpha)*(1-taue)*xb0/(r-mu0)-(1-taui)*c0/r)*(xi0/xb0)^beta2;
     
     phi0=1-(r*D0/((1-taui)*c0));
     
     Phi0=phi0;
     
   %%%%%Ͷ��ˮƽ%%%%%%%%%%%%
     y0=[1 0.8]';
     y0=fsolve(@(P)SB(P,Paras),y0);
 
     xi1=roundn(y0(1,1),-8);
     c1=roundn(y0(2,1),-8);
    
     xb1=(beta4/(beta4-1))*((r-mu1)/r)*c1;
     
     %%%%%%%������%%%%%%%%%%%%%%
     E1=(1-taue)*(xi1/(r-mu1)-c1/r+(c1/r-xb1/(r-mu1))*(xi1/xb1)^beta4);
     D1=(1-taui)*c1/r+((1-alpha)*(1-taue)*xb1/(r-mu1)-(1-taui)*c1/r)*(xi1/xb1)^beta4;
     
     phi1=1-(r*D1/((1-taui)*c1));
     
     Fg=(E1-phi1*(1-taui)*c1/r-(I-(1-taui)*c1/r))*(x/xi1)^beta3;
     
      %%%%%����Ͷ��ˮƽ%%%%%%%%%%%%
      y0=[1 0.8]';
      y0=fsolve(@(P)SEB0(P,Paras),y0);
 
     xis=roundn(y0(1,1),-8);  
      cs=roundn(y0(2,1),-8);  
    
     
          %%%%%%%������%%%%%%%%%%%%%%
       xbs=(beta4/(beta4-1))*((r-mu1)/r)*cs; 
     e1=(1-taue)*(xis/(r-mu1)-cs/r+(cs/r-xbs/(r-mu1))*(xis/xbs)^beta4);
     d1=(1-taui)*cs/r+((1-alpha)*(1-taue)*xbs/(r-mu1)-(1-taui)*cs/r)*(xis/xbs)^beta4;
     
     phis=1-(r*d1/((1-taui)*cs));
     Fs=(e1-phis*(1-taui)*cs/r-(I-(1-taui)*cs/r))*(x/xis)^beta3;
     
  %%%%%badtypepoolingͶ��ˮƽ%%%%%%%%%%%%
      y0=[1 1]';
      y0=fsolve(@(P)PB0(P,Paras),y0);
 
     xip0=roundn(y0(1,1),-8);  
    cp0=roundn(y0(2,1),-8); 
     
     
         
    
     
%                     %%%%%goodtypepoolingͶ��ˮƽ%%%%%%%%%%%%
%       y0=1;
%       y0=fsolve(@(P)PB1(P,Paras),y0);
%  
%      xip1=roundn(y0,-8);  
%      
     
     
                    %%%%%%%%%%%%%%�������poolingͶ��ˮƽ 
      xig=xi1:0.001:xi0;
   %  xig=xip0:0.001:xi0;
     n=length(xig);
     V=zeros(n,1);
     phip=zeros(n,1);
%      u=zeros(n,1);
     
     for j=1:n
      xip=xig(j);
      
      
       y0=1;
      y0=fsolve(@(P)Couponp(P,Paras,xip),y0);
  
    cpg=roundn(y0,-8); 
    
         
     %%%%%%%%%%%�Ʋ�ˮƽ%%%%%%%%%%%%%%
      xbp0=(beta2/(beta2-1))*((r-mu0)/r)*cpg;
     xbp1=(beta4/(beta4-1))*((r-mu1)/r)*cpg;
    
         %%%%%%%������%%%%%%%%%%%%%%
%      E0(j)=(1-taue)*(xip/(r-mu0)-cpg(j)/r+(cpg(j)/r-xbp0(j)/(r-mu0))*(xip0(j)/xb0)^beta2);
     D0(j)=(1-taui)*cpg/r+((1-alpha)*(1-taue)*xbp0/(r-mu0)-(1-taui)*cpg/r)*(xip/xbp0)^beta2;
     
%      E1(j)=(1-taue)*(xip/(r-mu1)-c/r+(c/r-xb1/(r-mu1))*(xip/xb1)^beta4);
     D1(j)=(1-taui)*cpg/r+((1-alpha)*(1-taue)*xbp1/(r-mu1)-(1-taui)*cpg/r)*(xip/xbp1)^beta4;
     
     phip(j)=1-(r*((1-pi0)*D0(j)+pi0*D1(j)))/((1-taui)*cpg);
    
    
    
    %%%%%%%��Ȩ��ֵ%%%%%%
    Em(j)=(1-taue)*(xip/(r-mu1)-cpg/r+(cpg/r-xbp1/(r-mu1))*(xip/xbp1)^beta4);
   % x=0.5; 
    
    V(j)=(Em(j)-phip(j)*(1-taui)*cpg/r-(I-(1-taui)*cpg/r))*(x/xip)^beta3;   
         
     u(j)= V(j);   
     end
     [L]=find(V==max(V));
  
     Xip=xig(L);
     Phip=phip(L);
     Fp=V(L);
%      u=u(L);
     
%       %%%%%%%high firm value%%%%%%%%%
% %        I=Paras(8);
%      Emh=(1-taue)*(Xip/(r-mu1)-c/r+(c/r-xb1/(r-mu1))*(Xip/xb1)^beta4);
%     Vh=(Emh-Phip*(1-taui)*c/r-(I-(1-taui)*c/r))*(x/Xip)^beta3;  
%      
%      
%      
%      %%%%%%%low firm value%%%%%%%%%
% %        I=Paras(8);
%      Eml=(1-taue)*(Xip/(r-mu0)-c/r+(c/r-xb0/(r-mu0))*(Xip/xb0)^beta2);
% %      M1=Phip*(1-taui)*c/r
% %      M2=(1-taui)*c/r
%     Vl=(Eml-Phip*(1-taui)*c/r-(I-(1-taui)*c/r))*(x/Xip)^beta1;  
%          
%      
     
     
     
     
     
%       %%%%%��СͶ��ˮƽ%%%%%%%%%%%%
%       y0=1;
%       y0=fsolve(@(P)SEB1(P,Paras),y0);
%  
%      xim=roundn(y0,-8);  
     
     
             
     
     
    %%%%%%%%%�ɱ�%%%%%%%%%%
    Cs=(Fg-Fs)/Fg;
    Cp=(Fg-Fp)/Fg;
   
     
     
     
     
     
     


end

function F=FB(P,Paras)
 xi0=P(1);
 c0=P(2);
%%%%%��������%%%%%%%%%%
  
    r=Paras(1);
    mu0=Paras(2);
    mu1=Paras(3);
    sigma=Paras(4);
    alpha=Paras(5);
    taue=Paras(6);
    taui=Paras(7);
    I=Paras(8);
    p=Paras(9);
    x=Paras(10);
    pi0=Paras(11);
    eta=Paras(12);
    
%     c=p*r/(1-taui);
    
   % psi=mu1-mu0;
    
    
    %%%%%%%%%%��%%%%%%%%
    beta1=1/2-mu0/sigma^2+sqrt((1/2-mu0/sigma^2)^2+2*r/sigma^2);
    beta2=1/2-mu0/sigma^2-sqrt((1/2-mu0/sigma^2)^2+2*r/sigma^2);
    
    xb0=(beta2/(beta2-1))*((r-mu0)/r)*c0;
%     p=(1-taui)*c/r+((1-alpha)*(1-taue)*xb0/(r-mu0)-(1-taui)*c/r)*(xi0/xb0)^beta2;

     F(1)=((taue-taui)*c0/r-I)*beta1+(beta1-1)*(1-taue)*xi0/(r-mu0)+(beta1-beta2)*((taui-taue)*c0/r-alpha*(1-taue)*xb0/(r-mu0))*(xi0/xb0)^beta2;
     
     F(2)=(1-taui)*c0/r+((1-alpha)*(1-taue)*xb0/(r-mu0)-(1-taui)*c0/r)*(xi0/xb0)^beta2-p;


end

function F=SB(P,Paras)
 xi1=P(1);
 c1=P(2);
%%%%%��������%%%%%%%%%%
%  global xb1 
    r=Paras(1);
    mu0=Paras(2);
    mu1=Paras(3);
    sigma=Paras(4);
    alpha=Paras(5);
    taue=Paras(6);
    taui=Paras(7);
    I=Paras(8);
    p=Paras(9);
    x=Paras(10);
    pi0=Paras(11);
    eta=Paras(12);
    
   % c=p*r/(1-taui);
    
   % psi=mu1-mu0;
    
    
    %%%%%%%%%%��%%%%%%%%
    beta3=1/2-mu1/sigma^2+sqrt((1/2-mu1/sigma^2)^2+2*r/sigma^2);
    beta4=1/2-mu1/sigma^2-sqrt((1/2-mu1/sigma^2)^2+2*r/sigma^2);
    
    xb1=(beta4/(beta4-1))*((r-mu1)/r)*c1;

     F(1)=((taue-taui)*c1/r-I)*beta3+(beta3-1)*(1-taue)*xi1/(r-mu1)+(beta3-beta4)*((taui-taue)*c1/r-alpha*(1-taue)*xb1/(r-mu1))*(xi1/xb1)^beta4;
  
     F(2)=(1-taui)*c1/r+((1-alpha)*(1-taue)*xb1/(r-mu1)-(1-taui)*c1/r)*(xi1/xb1)^beta4-p;

end

function F=SEB0(P,Paras)
 global  xb0  c0  Xi0 Phi0

xis=P(1);
 cs=P(2);
%%%%%��������%%%%%%%%%%

    r=Paras(1);
    mu0=Paras(2);
    mu1=Paras(3);
    sigma=Paras(4);
    alpha=Paras(5);
    taue=Paras(6);
    taui=Paras(7);
    I=Paras(8);
    p=Paras(9);
    x=Paras(10);
    pi0=Paras(11);
    eta=Paras(12);
    
    xi0=Xi0;
    phi0=Phi0;
    
%     c=p*r/(1-taui);
    
   % psi=mu1-mu0;
    
    
    %%%%%%%%%%��%%%%%%%%
    beta1=1/2-mu0/sigma^2+sqrt((1/2-mu0/sigma^2)^2+2*r/sigma^2);
    beta2=1/2-mu0/sigma^2-sqrt((1/2-mu0/sigma^2)^2+2*r/sigma^2);
    
       beta3=1/2-mu1/sigma^2+sqrt((1/2-mu1/sigma^2)^2+2*r/sigma^2);
    beta4=1/2-mu1/sigma^2-sqrt((1/2-mu1/sigma^2)^2+2*r/sigma^2);
    
  
    
    %%%%%%%��Ȩ��ֵ%%%%%%
      xbs=(beta2/(beta2-1))*((r-mu0)/r)*cs;
    Es=(1-taue)*(xis/(r-mu0)-cs/r+(cs/r-xbs/(r-mu0))*(xis/xbs)^beta2);
    
    xbg=(beta4/(beta4-1))*((r-mu1)/r)*cs;
    Ds=(1-taui)*cs/r+((1-alpha)*(1-taue)*xbg/(r-mu1)-(1-taui)*cs/r)*(xis/xbg)^beta4;
    phis=1-Ds/((1-taui)*cs/r);
    
    E0=(1-taue)*(xi0/(r-mu0)-c0/r+(c0/r-xb0/(r-mu0))*(xi0/xb0)^beta2);
    
    
    

   F(1)=Es-phis*(1-taui)*cs/r-(I-(1-taui)*cs/r)-(E0-phi0*(1-taui)*c0/r-(I-(1-taui)*c0/r))*(xis/xi0)^beta1;
   
   F(2)=(1-phis)*(1-taui)*cs/r-p;


end

% function F=SEB1(P,Paras)
%  xim=P;
% %%%%%��������%%%%%%%%%%
%  global xb0 xb1 phi0 phi1 Phip xi0 xi1  Xip
%     r=Paras(1);
%     mu0=Paras(2);
%     mu1=Paras(3);
%     sigma=Paras(4);
%     alpha=Paras(5);
%     taue=Paras(6);
%     taui=Paras(7);
%     I=Paras(8);
%     p=Paras(9);
%     x=Paras(10);
%     pi0=Paras(11);
%     eta=Paras(12);
%     
%     c=p*r/(1-taui);
%     
%    % psi=mu1-mu0;
%     
%     
%     %%%%%%%%%%��%%%%%%%%
%     beta3=1/2-mu1/sigma^2+sqrt((1/2-mu1/sigma^2)^2+2*r/sigma^2);
%     beta4=1/2-mu1/sigma^2-sqrt((1/2-mu1/sigma^2)^2+2*r/sigma^2);
%     
%     %%%%%%%��Ȩ��ֵ%%%%%%
%     Em=(1-taue)*(xim/(r-mu1)-c/r+(c/r-xb1/(r-mu1))*(xim/xb1)^beta4);
%     E1=(1-taue)*(Xip/(r-mu1)-c/r+(c/r-xb1/(r-mu1))*(Xip/xb1)^beta4);
%    % x=0.5; 
%     
%     
%     
% 
%    F=Em-phi1*(1-taui)*c/r-(I-(1-taui)*c/r)-(E1-Phip*(1-taui)*c/r-(I-(1-taui)*c/r))*(xim/xi0)^beta3;
% 
% 
% end

function F=PB0(P,Paras)
 global  xb0  c0  Xi0 Phi0


xip0=P(1);
 cp0=P(2);
%%%%%��������%%%%%%%%%%
    r=Paras(1);
    mu0=Paras(2);
    mu1=Paras(3);
    sigma=Paras(4);
    alpha=Paras(5);
    taue=Paras(6);
    taui=Paras(7);
    I=Paras(8);
    p=Paras(9);
    x=Paras(10);
    pi0=Paras(11);
    eta=Paras(12);
    
%     c=p*r/(1-taui);
    
   % psi=mu1-mu0;
    xi0=Xi0;
    phi0=Phi0;
    
    %%%%%%%%%%��%%%%%%%%
    beta1=1/2-mu0/sigma^2+sqrt((1/2-mu0/sigma^2)^2+2*r/sigma^2);
    beta2=1/2-mu0/sigma^2-sqrt((1/2-mu0/sigma^2)^2+2*r/sigma^2);
       
    beta3=1/2-mu1/sigma^2+sqrt((1/2-mu1/sigma^2)^2+2*r/sigma^2);
    beta4=1/2-mu1/sigma^2-sqrt((1/2-mu1/sigma^2)^2+2*r/sigma^2);
    
    %%%%%%%%%%%�Ʋ�ˮƽ%%%%%%%%%%%%%%
    %%%%�Ʋ�ˮƽ%%%%%
    xbp0=(beta2/(beta2-1))*((r-mu0)/r)*cp0;
    xbp1=(beta4/(beta4-1))*((r-mu1)/r)*cp0;
    
         %%%%%%%������%%%%%%%%%%%%%%
%      E0=(1-taue)*(xip0/(r-mu0)-cp0/r+(cp0/r-xbp0/(r-mu0))*(xip0/xb0)^beta2);
     D0=(1-taui)*cp0/r+((1-alpha)*(1-taue)*xbp0/(r-mu0)-(1-taui)*cp0/r)*(xip0/xbp0)^beta2;
     
%      E1=(1-taue)*(xip0/(r-mu1)-c/r+(c/r-xb1/(r-mu1))*(xip0/xb1)^beta4);
     D1=(1-taui)*cp0/r+((1-alpha)*(1-taue)*xbp1/(r-mu1)-(1-taui)*cp0/r)*(xip0/xbp1)^beta4;
     
     phip0=1-(r*((1-pi0)*D0+pi0*D1))/((1-taui)*cp0);
    
    
    
    %%%%%%%��Ȩ��ֵ%%%%%%
    Es=(1-taue)*(xip0/(r-mu0)-cp0/r+(cp0/r-xbp0/(r-mu0))*(xip0/xbp0)^beta2);
    E0=(1-taue)*(xi0/(r-mu0)-c0/r+(c0/r-xb0/(r-mu0))*(xi0/xb0)^beta2);
    

   F(1)=Es-phip0*(1-taui)*cp0/r-(I-(1-taui)*cp0/r)-(E0-phi0*(1-taui)*c0/r-(I-(1-taui)*c0/r))*(xip0/xi0)^beta1;
   
   F(2)=(1-phip0)*(1-taui)*cp0/r-p;


end

function F=Couponp(P,Paras,xip)

 cpg=P;
%%%%%��������%%%%%%%%%%
%  global xb0 xb1 phi0 phi1 xi0 xi1 phip0 c0
    r=Paras(1);
    mu0=Paras(2);
    mu1=Paras(3);
    sigma=Paras(4);
    alpha=Paras(5);
    taue=Paras(6);
    taui=Paras(7);
    I=Paras(8);
    p=Paras(9);
    x=Paras(10);
    pi0=Paras(11);
    eta=Paras(12);
    
    xi=xip;
%     c=p*r/(1-taui);
    
   % psi=mu1-mu0;
    
    
    %%%%%%%%%%��%%%%%%%%
    beta1=1/2-mu0/sigma^2+sqrt((1/2-mu0/sigma^2)^2+2*r/sigma^2);
    beta2=1/2-mu0/sigma^2-sqrt((1/2-mu0/sigma^2)^2+2*r/sigma^2);
       
    beta3=1/2-mu1/sigma^2+sqrt((1/2-mu1/sigma^2)^2+2*r/sigma^2);
    beta4=1/2-mu1/sigma^2-sqrt((1/2-mu1/sigma^2)^2+2*r/sigma^2);
    
    %%%%%%%%%%%�Ʋ�ˮƽ%%%%%%%%%%%%%%
    %%%%�Ʋ�ˮƽ%%%%%
    xbp0=(beta2/(beta2-1))*((r-mu0)/r)*cpg;
    xbp1=(beta4/(beta4-1))*((r-mu1)/r)*cpg;
    
         %%%%%%%������%%%%%%%%%%%%%%
%      E0=(1-taue)*(xip0/(r-mu0)-cp0/r+(cp0/r-xbp0/(r-mu0))*(xip0/xb0)^beta2);
     D0=(1-taui)*cpg/r+((1-alpha)*(1-taue)*xbp0/(r-mu0)-(1-taui)*cpg/r)*(xi/xbp0)^beta2;
     
%      E1=(1-taue)*(xip0/(r-mu1)-c/r+(c/r-xb1/(r-mu1))*(xip0/xb1)^beta4);
     D1=(1-taui)*cpg/r+((1-alpha)*(1-taue)*xbp1/(r-mu1)-(1-taui)*cpg/r)*(xi/xbp1)^beta4;
     
     phip1=1-(r*((1-pi0)*D0+pi0*D1))/((1-taui)*cpg);
    
    
    
%     %%%%%%%��Ȩ��ֵ%%%%%%
%     Es=(1-taue)*(xip0/(r-mu0)-cp0/r+(cp0/r-xbp0/(r-mu0))*(xip0/xbp0)^beta2);
%     E0=(1-taue)*(xi0/(r-mu0)-c0/r+(c/r-xb0/(r-mu0))*(xi0/xb0)^beta2);
%     
% 
%    F(1)=Es-phip0*(1-taui)*cp0/r-(I-(1-taui)*cp0/r)-(E0-phi0*(1-taui)*c0/r-(I-(1-taui)*c0/r))*(xip0/xi0)^beta1;
   
   F(2)=(1-phip1)*(1-taui)*cpg/r-p;


end

% function F=PB1(P,Paras)
%  xip1=P;
% %%%%%��������%%%%%%%%%%
%  global xb0 xb1 phi0 phi1 phis xi0 xi1 xis xim phip1
%     r=Paras(1);
%     mu0=Paras(2);
%     mu1=Paras(3);
%     sigma=Paras(4);
%     alpha=Paras(5);
%     taue=Paras(6);
%     taui=Paras(7);
%     I=Paras(8);
%     p=Paras(9);
%     x=Paras(10);
%     pi0=Paras(11);
%     eta=Paras(12);
%     
%     c=p*r/(1-taui);
%     
%    % psi=mu1-mu0;
%     
%     
%     %%%%%%%%%%��%%%%%%%%
%     beta1=1/2-mu0/sigma^2+sqrt((1/2-mu0/sigma^2)^2+2*r/sigma^2);
%     beta2=1/2-mu0/sigma^2-sqrt((1/2-mu0/sigma^2)^2+2*r/sigma^2);
%        
%     beta3=1/2-mu1/sigma^2+sqrt((1/2-mu1/sigma^2)^2+2*r/sigma^2);
%     beta4=1/2-mu1/sigma^2-sqrt((1/2-mu1/sigma^2)^2+2*r/sigma^2);
%     
%         %%%%%%%%%%%�Ʋ�ˮƽ%%%%%%%%%%%%%%
%      %%%%�Ʋ�ˮƽ%%%%%
%     xb0=(beta2/(beta2-1))*((r-mu0)/r)*c;
%     xb1=(beta4/(beta4-1))*((r-mu1)/r)*c;
%     
%          %%%%%%%������%%%%%%%%%%%%%%
%      E0=(1-taue)*(xip1/(r-mu0)-c/r+(c/r-xb0/(r-mu0))*(xip1/xb0)^beta2);
%      D0=(1-taui)*c/r+((1-alpha)*(1-taue)*xb0/(r-mu0)-(1-taui)*c/r)*(xip1/xb0)^beta2;
%      
%      E1=(1-taue)*(xip1/(r-mu1)-c/r+(c/r-xb1/(r-mu1))*(xip1/xb1)^beta4);
%      D1=(1-taui)*c/r+((1-alpha)*(1-taue)*xb1/(r-mu1)-(1-taui)*c/r)*(xip1/xb1)^beta4;
%      
%      phip1=1-(r*((1-pi0)*D0+pi0*D1))/((1-taui)*c);
%     
%     
%     
%     %%%%%%%��Ȩ��ֵ%%%%%%
%     Em=(1-taue)*(xip1/(r-mu1)-c/r+(c/r-xb1/(r-mu1))*(xip1/xb1)^beta4);
%     E1=(1-taue)*(xis/(r-mu1)-c/r+(c/r-xb1/(r-mu1))*(xis/xb1)^beta4);
%     
% 
%    F=(Em-phip1*(1-taui)*c/r-(I-(1-taui)*c/r))*(1/xip1)^beta3-(E1-phis*(1-taui)*c/r-(I-(1-taui)*c/r))*(1/xis)^beta3;
% 
% end
% 
% 
% 
